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Abstract. This paper studies the stability properties of a two dimensional rela¬ 
tive velocity scheme for the Navier-Stokes equations. This scheme inspired by the 
cascaded scheme has the particularity to relax in a frame moving with a velocity 
field function of space and time. Its stability is studied first in a linear context 
then on the non linear test case of the Kelvin-Helmholtz instability. The link with 
the choice of the moments is put in evidence. The set of moments of the cascaded 
scheme improves the stability of the d’Humieres scheme for small viscosities. On 
the contrary, a relative velocity scheme with the usual set of moments deteriorates 
the stability. 


Introduction 

The lattice Boltzmann schemes have been successfully used for the simulation of the 
compressible Navier-Stokes equations in two or three dimensions [?, ?, ?, ?]. This 
method aims to mimic the microscopic behaviour in order to simulate some macro¬ 
scopic problems. The algorithm consists in evaluating some particle distributions. 
The particles, moving from node to node of a lattice, undergo a phase of collision 
and a phase of transport. Different collision operators have been proposed for the 
simulation of the Navier-Stokes equations. The simplest one is the single relaxation 
time operator [?, ?, ?, ?] also called BGK. An alternative called the multiple relax¬ 
ation times (MRT) operator [?, ?] has been proposed. During the collision, some 
moments, linear combinations of the particle distributions, relax towards the equi¬ 
librium with a priori different velocities. It contains the particularity to offer more 
degrees of freedom to fix the different parameters as the viscosities. The multiple 
relaxation times approach is thus more flexible than the BGK. Both schemes have 
been well studied particularly in terms of stability [?]. They still encounter some 
instability features as the viscosities tend to zero that limits high Reynolds number 
simulations. 

In 2006, a cascaded scheme improving the stability for low viscosities has been pre¬ 
sented [?]. Its relaxation occurs in a frame moving with the fluid velocity. To un¬ 
derstand the positive features of this scheme, a general notion of relative velocity 
schemes was defined [?]. Their relaxation is made for a set of moments depending on 
a velocity held function of space and time that is the velocity huid for the cascaded 
scheme [?] and zero for the d’Humieres scheme [?]. These relative velocity schemes 
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are not restricted to the simulation of the Navier-Stokes equations: they are defined 
for an arbitrary number of conservation laws. Their consistency has already been 
studied for one and two conservation laws [?, ?, ?] but the same reasoning holds for 
an arbitrary number of conservation laws. 

The purpose of this contribution is to present some numerical stability results of the 
two dimensional nine velocities (D 2 Q 9 ) relative velocity scheme for the compressible 
Navier-Stokes equations. We want to characterize the influence of the relative velocity 
and the link with the moments choice: the polynomials defining the moments of the 
cascaded scheme are different from the usual ones and may act on the stability. In 
a first part, we recall the basis of the relative velocity schemes. We then present the 
relative velocity D 2 Q 9 we are interested in. The second part exhibits the results of 
stability, first in a linear context (L^ von Neumann notion) and then for a non linear 
test case, the Kelvin-Helmholtz instability. It puts in evidence the link between the 
relative velocity, the choice of the polynomials defining the moments and the stability. 

1. Description of the scheme 

We first introduce the relative velocity scheme for an arbitrary number of dimensions 
and velocities. We then particularize it to the case of two dimensions and nine 
velocities. 

1.1. The relative velocity scheme 

This section presents the derivation of the relative velocity lattice Boltzmann schemes 
introduced in [?] and inspired by the cascaded scheme [?]. Let £ be a cartesian lat¬ 
tice in d dimensions with a typical mesh size Ax. The time step At is linked to the 
space step by the acoustic scaling At = Ax/A for A G M the velocity scale. We in¬ 
troduce V = (r)g,..., a set of q velocities of M'^. This defines the scheme called 
DrfQg. We assume that for each node x of the lattice £, and each v- in V, the point 
X + VjAt is still a node of C. The D^Qg scheme computes a particle distribution 
f = {/q, ..., fg_i) on the lattice C at discrete values of time. An iteration of the 
scheme consists in two phases: the relaxation that is non linear and local in space, 
and the linear transport solved exactly by a characteristic method. 

The relaxation phase reads more easily in a moments basis using the d’Humieres 
framework [?]. A velocity field u{x,t) that depends on space and time being given, 
we define the matrix of moments M{u) by 

Mi^)kj = ^ki^j - 0 ^ k,j ^ q-1, 

where {Pq, ..., Pg_i) are some polynomials of ..., X^]. This matrix of moments, 

supposed to be invertible, defines the moments m{u) = (mg(?t),..., mg_^{u)) by the 
relation 

(1) m{u) = M{u)f, 

where m^{u) is the moment. 



ON THE STABILITY OF A RELATIVE VELOCITY LATTICE BOLTZMANN SCHEME 


3 


The relative velocity schemes use a diagonal relaxation phase in the shifted moments 
basis 

(2) ml{u) = 0 ^ A: ^ q-l, 

where m^{u) is the A:**' moment at equilibrium and the relaxation parameter as¬ 
sociated with the moment for 0 ^ A: ^ <?~1- Some of these moments are conserved 
by the relaxation: they are associated with relaxation parameters equal to zero. The 
equilibrium derives from a vector of distribution functions at the equilibrium 
independent of the velocity field u, only dependent of the conserved moments. 

(3) m^\u) = M{u)f^^. 

The inverse of the matrix of moments is used to return to the distributions 

(4) /* = M~^{u)m*{u). 

The transport phase spreads the particle distributions on the neighbouring nodes 

fj{x,t + At) = f*{x- VjAt,t), q-l. 

This framework embeddes the d’Humieres scheme for u equal to 0 and the cascaded 
scheme for u equal to the fluid velocity and a particular set of moments [?]. In the 
following, when u is specified, we call the associated relative velocity scheme the 
scheme relative to u. 

1.2. The study framework: the relative velocity D 2 Q 9 scheme 

The purpose of this section is to introduce the scheme whose stability properties are 
investigated: the relative velocity D 2 Q 9 scheme with two conservation laws on the 
density and the momentum 

P = '^ fj^ 1'^ = Y1 l^a^d. 

3 3 

for the compressible Navier-Stokes equations. We expose its features and the differ¬ 
ent degrees of freedom used to check its stability. We put a particular attention on 
the definition of the moments. 

For this two-dimensional scheme, nine velocities are involved: they are defined by 
r; = {(0,0), (A, 0), (0, A), (-A, 0), (0, -A), (A, A), (-A, A), (-A, -A), (A, -A)}, 
with A G M the velocity scale. These velocities are also represented on the figure 

We need to deal with the set of the moments and the equilibrium to completely 
characterize the scheme. Historically the D 2 Q 9 scheme has been mainly used with 
the following set of moments 

(5) 1, X, Y, - Y'^,XY, X(X^ + Y^), Y(X^ + Y^), (X^ + Y^f, 

or its orthogonalized analogue for the simulation of the Navier-Stokes equations [?]. 
They have been chosen because of their physical meaning: they involve the density, 
the momentum, the energy, the diagonal and off-diagonal components of the stress 
tensor, the heat flux and the square of the energy. Nevertheless, the D 2 Q 9 cascaded 


4 


FRANgOIS DUBOIS, TONY FEVRIER, AND BENJAMIN GRAILLE 


A 



Figure 1 . The D 2 Q 9 velocities. 


scheme [?], that seems to improve the stability at low viscosities, has brought to light 
an other set of moments given by 

(6) i,x, Y, ^2 _ y^,xy, xy^, yx"^, 

This scheme has been written as a relative velocity scheme for the moments (|^ and 
the fluid velocity [?]. The relaxation of ([^ and Q are equivalent in the d’Humieres 
framework corresponding to it = 0 (section |2.2[ ). However, this is not true any more 
when u is different from 0. This point naturally leads to the question: does the set of 
moments have an influence on the stability properties of the relative velocity scheme? 
Giving some experimental rudiments of an answer is the purpose of this study. 

That’s why we introduce two sets of moments tuned by a parameter a G M: 

(7) 1, y, y, x^+Y^, x^-Y^, XY, x{ax‘^+Y^),Y{x'^+aY‘^), |(x^+y^)+x2y2, 

and 

(8) 1 , y, y, + y^, - y^ xy, xy^ + a{x‘^ + y^), yx^ + a(x2 + y^), x^y^. 

The moments ([^ generalize those of the D2Q9 cascaded scheme corresponding to 
a = 0 [?] given by and the ones associated with a = 1 defined by ©■ The intro¬ 
duction of a results from the will not to restrict the study to two sets of moments. 
This allows also to understand the impact on the stability of the X^ component when 
a moves from 0 to 1. The choice of the moments ([^, even if it seems strange because 
it mixes some second and third order polynomials, improves the understanding of 
the differences of stability between ([^ and (|^ for the relative velocity D2Q9 scheme. 
Taking a = 0 also recovers the cascaded moments. 
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The equilibrium may also have an influence on the stability. That’s why, denoting 
u = q!p the fluid velocity, we introduce 


( 9 ) 

and 


u) = pujj 1 H-^ + 


u.v- 




2Cn 


\u\ 


2Cn 


0 ^ ^ 8 , 


( 10 ) 


fPip,u) = pUJj 


1 + 


u.v- 


+ 


(U.v- 


2Cn 


\uy 


2c" 


+ 


{u.V-f 

n 6 

6Cn 


\u\ 


’{u.v-) dJu^)‘^{uy)‘^ 


d. = 1/2, j = 1,...,4, 


2co 

dj = -1’ 


+ 


0 ^ ^ 8 , 


j = 5 ,..., 8, respectively corre- 


where dg = “ 1 / 4 , 

spending to the second order truncated equilibrium [?] and to the product equilibrium 
used for the D2Q9 cascaded scheme [?]. The product equilibrium corresponds to the 
fourth order truncation of the maxwellian equilibrium. Both equilibria allow to sim¬ 
ulate the compressible Navier-Stokes equations whatever the velocity field u. Indeed, 
the second order equivalent equations of the relative velocity schemes are indepen¬ 
dent of u [?] and the Navier-Stokes equations are recovered by the D2Q9 d’Humieres 
scheme {u = 0) at the second order for small Mach numbers [?]. Let’s note that the 
simulations of the section have been also made for the incompressible analogue of 
(|^ and (10) used in [?]: same trends as those presented in the section]^ are obtained. 


We choose to work with several two relaxation times schemes (TRT) to understand 
the role of each polynomial of the moments: the one given by 

(f f) ^ {^ei 1 ^ 1 ^ei ^ei ^e) 1 

called TRTi and the one given by 

(12) s (Sg, Sg, Sg, Sg), 

called TRT2 where Sg, Sp G M. Note that the TRTi and the TRT2 differ from the 
TRT schemes defined in [?] and based on the symmetry of the lattice. If all the 
relaxation parameters are identical, we recover the BGK scheme [?]. 


Four degrees of freedom are tunable in this section: the moments, the vector of the 
relaxation parameters s, the velocity field u and the equilibrium. The link between 
these parameters and their influence on the stability is studied in the following. 


2. Experimental study of linear stability 


In this section, we study the linear von Neumann L?' stability of the relative velocity 
D2Q9 scheme defined in the section 1.2 The influence of the moments according to 


the choice of the velocity field u is the keypoint of the section. Our first interest goes 
to the moments (|^ and (|^, respectively corresponding to a = 0 and a = 1 in ([^, 
because they are usually chosen by the community [?, ?, ?, ?]. The first subsection 
compares those two sets according to the velocity field parameter. We show that 
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taking u = u, the velocity of the fluid, improves the stability if the moments Q 
are chosen, deteriorates it if the set (|^ is taken. The second subsection answers the 
following question: what is the better choice of moments (of a) in terms of stability? 
A range of a is proposed and a = 0 is showed to be the most stable choice. 


2.1. Methodology: the von Neumann stability 

The study of the relative velocity D 2 Q 9 scheme is based on the von Neumann 
stability. This notion being adapted to linear contexts, we linearize the equilibria ([^ 


and (10) around a velocity V = £ M^, 0 G M. Thus there exists a matrix E 


so that 

= Ef. 

Using (HHH the linearized relaxation phase of the relative velocity schemes reads 
r = (/ + M{u)-^DM{u){E - I))f, 


where D = diag(s) is the diagonal matrix of the relaxation parameters. This expres¬ 
sion holds for each node x of the lattice, the relaxation being local in space. One can 
deduce the expression of the distribution after an iteration thanks to the transport 
phase 

f-{x,t + At) = [{I + M{u)-^DM{u){E-I))f]j{x-v.At,t), x e C, teR. 

In the Fourier space, the transport operator becomes local in space and is represented 
by the diagonal matrix A whose diagonal components are given by , 0 ^ j ^ 

8 . We can then define the amplification matrix L(u) = L{u,V,k ,s,a) = A{I 
M{u)~^DM{u){E — /)), for k,V,u £ M^, s G M®, a G M, G characterizing 
a time iteration of the scheme in the Fourier space 

f{k,t + At) = L{u)f{k,t), t G M, 

where f is the Fourier transform of f. We want to determine the quantity 


(13) max{|'U|, max r(L{u)) ^ 1}, 

fceiR2 

for some parameters s, u, a, a direction of linearization 0 G M and r{L{u)) the 
spectral radius of L[u). It characterizes the set of the linearization velocities V for 
which the scheme verifies the necessary condition of stability max r(L(u)) ^ 1. 

feeK2 


2.2. Comparison between the d’Humieres scheme and the scheme relative 
to the linearization velocity 

We show that the schemes relative to u = V can improve or deteriorate the lin¬ 
ear stability compared to the d’Humieres scheme. The stability behaviour depends 
strongly on the choice of the moments. 


We compare the schemes relative to S = 0 and u = V for the two sets of moments 
(|^ and ([^: we have a = 0,1 in ([^. We here restrict to the second order truncated 
equilibrium linearized around V. The variable of comparison is the largest stable 
velocity V (13) for a linearization direction 0 equal to 0. We choose to deal with 
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n 

m 0 

1 

2 

3 

4 

5 

6 

7 

0 

0.42 

0.41 

0.34 

0.26 

0.20 

0.15 

0.11 

0.08 

1 

0.42 

0.41 

0.36 

0.30 

0.23 

0.18 

0.13 

0.09 

2 

0.31 

0.34 

0.34 

0.32 

0.28 

0.23 

0.17 

0.13 

3 

0.21 

0.28 

0.32 

0.30 

0.25 

0.22 

0.18 

0.15 

4 

0.14 

0.21 

0.28 

0.26 

0.22 

0.18 

0.16 

0.13 

5 

0.10 

0.16 

0.22 

0.23 

0.20 

0.17 

0.13 

0.11 

6 

0.07 

0.12 

0.17 

0.20 

0.18 

0.16 

0.12 

0.11 

7 

0.05 

0.08 

0.12 

0.17 

0.16 

0.15 

0.11 

0.11 


Table 1. Highest stable V - 

= (W 

',0) in A units for the d’Humieres 


scheme (7t = 0), with a = 0 oi 

■ a = 

1 , of equilibrium 0. Sg 

= 2-2“™ 


and = 2 — 2 









n 

ra 

0 

1 

2 

3 

4 

5 

6 

7 

0 


0.42 

0.42 

0.40 

0.37 

0.33 

0.28 

0.24 

0.21 

1 


0.42 

0.41 

0.37 

0.34 

0.33 

0.30 

0.27 

0.23 

2 


0.42 

0.36 

0.34 

0.33 

0.29 

0.25 

0.22 

0.19 

3 


0.36 

0.34 

0.33 

0.30 

0.25 

0.21 

0.18 

0.16 

4 


0.32 

0.32 

0.30 

0.26 

0.22 

0.18 

0.16 

0.13 

5 


0.29 

0.30 

0.27 

0.24 

0.21 

0.17 

0.13 

0.11 

6 


0.26 

0.28 

0.24 

0.21 

0.18 

0.16 

0.12 

0.11 

7 


0.23 

0.26 

0.22 

0.19 

0.16 

0.15 

0.11 

0.11 


Table 2. Highest stable V = ( V ^, 0 ) in A units for the scheme rel¬ 
ative to u = V with a = 0 of equilibrium = 2 — 2“™" and 

s, = 2 - 2-T 


the TRTi (11) for Sg = 2 — 2“"* and Sj, = 2 — 2“"' where m, n, G N, 0 ^ m, n ^ 7. 
The parameters Sg and respectively tune the bulk and the shear viscosities of the 
Navier-Stokes equations. This choice of parameters allows to study the zero viscos¬ 
ity limit by increasing m or/and n. The table deals with the d’Humieres scheme 
for both sets of moments, the values for those two sets being identical. The tables 
and 1^ give analogous results for the scheme relative to u = V: they correspond 
respectively to the moments with a = 0 ([^ and a = 1 (|^. 


We notice the importance of the choice of the moments for the schemes relative to 
u = V: stability areas are the biggest for a = 0 0 (table 0 and the smallest for 
a = 1 0 (table 0 whatever the choice of s. The d’Humieres scheme (S = 0, table 
0 has smaller stability areas than the scheme relative to u = V with a = 0 and 
bigger than the one with a = 1. 
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n 

ra 

0 

1 

2 

3 

4 

5 

6 

7 

0 


0.42 

0.42 

0.27 

0.18 

0.12 

0.08 

0.06 

0.04 

1 


0.42 

0.41 

0.37 

0.31 

0.20 

0.13 

0.09 

0.06 

2 


0.24 

0.36 

0.34 

0.32 

0.27 

0.21 

0.14 

0.09 

3 


0.15 

0.29 

0.33 

0.30 

0.25 

0.21 

0.18 

0.14 

4 


0.10 

0.19 

0.29 

0.26 

0.22 

0.18 

0.16 

0.13 

5 


0.07 

0.12 

0.20 

0.29 

0.21 

0.17 

0.13 

0.11 

6 


0.05 

0.09 

0.14 

0.20 

0.18 

0.16 

0.12 

0.11 

7 


0.03 

0.06 

0.09 

0.14 

0.16 

0.15 

0.11 

0.11 


Table 3. Highest stable V = (V^,0) in A units for the scheme rel¬ 
ative to u = V with a = 1 of equilibrium ([^. Sg = 2 — 2“™" and 
s, = 2 - 2-T 


The scheme relative to u = V with a = 0 provides the most important gain com¬ 
pared to the d’Humieres scheme when Sg or is close to 2 and the other is far from 2 
(for one small and one large viscosity). Instead these areas are the most deteriorated 
when a is equal to 1. When Sg and are close, the scheme presents stability areas 
nearly independent of u. The case Sg = corresponds to the BGK scheme: the 
velocity field u does not play any role since ^ is verified. 

Finally, for the d’Humieres scheme, the results are independent of the choice of the 
moments. Indeed, the third order moment X{X'^ is equivalent to X^X + XY'^ 

on the velocity network [?]: its relaxation is then equivalent to the relaxation of XY^, 
X being a conserved component. The same reasoning holds for the symmetrical mo¬ 
ment. The moment {X^ +Y^)‘^ is equal to X‘^Y^ + + Y"^) on the velocity 

network. Its relaxation is equivalent to relax X Y because X + Y and X Y are 
both in the eigenspace related to Sg for the TRTi. 


All these trends are independent of the direction of linearization 9: the results are 
similar for 6 = tt/S, vr/d, vr/S. 


We now do the same job for the product equilibrium (10) restricting the study to 
the choice a = 0. Note that the combination of these moments and this equilibrium 
corresponds to the D2Q9 cascaded scheme [?, ?]. The TRTi (11), tuned by Sg and 
s^, and the direction 9 = 0 are chosen. The table is about the d’Humieres scheme, 
the table corresponds to the scheme relative to u = V. 


The results are analogous to the ones associated with the equilibrium ([^ . The scheme 
relative to u = V has bigger linear stability areas than the d’Humieres scheme 
{u = 0) when a = 0. The gain is more important when the relaxation parameters 
are far from each other. The velocity field impact is lightened when Sg and are close. 
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n 

m 0 

1 

2 

3 

4 


5 

6 

7 

0 

0.42 

0.42 

0.39 

0.32 

0.24 

0.16 

0.11 

0.07 

1 

0.42 

0.42 

0.41 

0.38 

0.31 

0.20 

0.14 

0.09 

2 

0.42 

0.42 

0.41 

0.40 

0.38 

0.30 

0.20 

0.14 

3 

0.26 

0.41 

0.40 

0.39 

0.37 

0.32 

0.28 

0.20 

4 

0.16 

0.28 

0.38 

0.36 

0.33 

0.29 

0.24 

0.21 

5 

0.10 

0.18 

0.29 

0.33 

0.31 

0.26 

0.21 

0.19 

6 

0.07 

0.12 

0.19 

0.30 

0.29 

0.24 

0.20 

0.18 

7 

0.05 

0.08 

0.13 

0.20 

0.28 

0.23 

0.19 

0.17 


Table 4. Highest stable V 

= (^", 

0) in A units for the d’Humieres 


scheme (u = 0), with a = 0 

of equilibrium ( 

10 

). Sg 

= 2- 

2-”^ and 


A = 2 - 2-T 










n 

m 

0 

1 

2 

3 

4 

5 

6 

7 

0 


0.42 

0.42 

0.42 

0.42 

0.35 

0.30 

0.26 

0.23 

1 


0.42 

0.42 

0.42 

0.41 

0.39 

0.35 

0.32 

0.28 

2 


0.42 

0.41 

0.41 

0.40 

0.40 

0.35 

0.31 

0.29 

3 


0.41 

0.41 

0.40 

0.39 

0.36 

0.30 

0.27 

0.23 

4 


0.40 

0.40 

0.39 

0.36 

0.33 

0.28 

0.23 

0.20 

5 


0.35 

0.37 

0.36 

0.33 

0.31 

0.26 

0.21 

0.18 

6 


0.31 

0.33 

0.33 

0.31 

0.29 

0.25 

0.20 

0.17 

7 


0.28 

0.30 

0.31 

0.29 

0.27 

0.24 

0.19 

0.17 


Table 5. Highest stable V = {V^,0) in A units for the scheme rel¬ 
ative lo u = V with a = 0 of equilibrium (10). Sg = 2 — 2 “™ and 
&, = 2 - 2 -". 


We finally assess the influence of the equilibrium on the linear stability. Whatever 
the choice of u and s, the equilibrium (10) provides bigger stability areas than the 
truncated equilibrium ([^. Particularly, the BGK scheme associated with (10) is 
more stable than the one corresponding to (§. 


As a conclusion, the most important fact of the study is the following: the scheme 
relative to u = V for a = 0 is more stable than for a = 1. Instead choosing a scheme 
relative tou = V with a “inappropriate” set of moments can deteriorate the stability. 


2.3. Influence of the choice of the moments on the stability 

The previous section has studied the stability of the relative velocity and d’Humieres 
schemes for two choices of a. This parameter seems to be crucial for the relative 
velocity schemes. The purpose of this section is to see more precisely its influence 
on the stability. It studies the stability properties of the schemes relative to S = 0 
and u = V for a bigger range of a. We show numerically and justify that a = 0 
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-0.8 -0.6 -0.4 -0.2 


m=0,n=3 

m=0,n=7 

m=3,n=0 

m=7,n=0 

m=7,n=7 


0.2 0.4 0.6 0.8 


Figure 2. Draw of as a function of a for the d’Humieres scheme 
with the moments Q. Left: TRTi, Sg = 2 — 2“™' and = 2 — 2“”. 
Right: TRTa, Sg = 2 - 2"”^ and = 2 - 2"". 


constitutes the better choice of moments. 


We are interested in the stability of the relative velocity schemes for both sets of 
moments ([^ and (|^: the discussion carries on the choice of the parameter a G M 
characterizing these moments. Both equilibria leading to the same trends, we focus 
on the truncated one (|^ linearized around V = (R*, 0) G M^. Two sets of relaxation 
parameters s are used: the TRTi (11) and the TRT2 (12) where Sg = 2 — 2“”^ and 
= Sp = 2 — 2“"' with {m, n) = (0, 3), (3, 0), (0, 7), (7,0), (7, 7). The quantity (13) is 
drawn as a function of a in [—1,1]. A negative value of (13) means that the scheme 
is unstable for all V. 


We first focus on the d’Humieres scheme corresponding to = 0. The figures and 
[^represent respectively the draws associated with the moments (|^ and (|^. On each 
figure, the left draw is associated with the TRTi and the right one with the TRT2. 

For the moments (Q, the draws are independent of a whatever the TRT chosen and 
s. For the moments ([^, the draw corresponding to the TRTi is independent of a 
unlike the TRT2. The figure associated with the TRT2 induces to choose a = 0: 
it corresponds to the maximum of the curve and the stability area decreases as |a| 
increases. As expected, the draw for m = n = 7 corresponding to a BGK scheme is 
constant in a. We notice that a = 0 belongs to the set of a maximizing the stability 
whatever the draw. 
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Figure 3. Draw of as a function of a for the d’Humieres scheme 
with the moments Left: TRTi, Sg = 2 — 2“™' and = 2 — 2“”. 

Right: TRTa, Sg = 2 - 2"”^ and = 2 - 2-". 

We can exhibit the origin of the dependence or independence on a. Let’s consider 
the moments ([^. For the d’Humieres scheme, the relaxation of these moments is 
independent of a. The last three moments of ([^ are 

+ XY^, aY^ + + Y^) + X^Y^^. 

Knowing that X^ = X'^X on the velocity set [?], the scheme is unchanged if we 
replace them by 

+ XY^, X^aY + X^Y, ^{X^ + Y^) + X^Y^. 

Relaxing the moments Q is then equivalent to relax the same moments for a = 0. 
Indeed, X and Y are associated with some conserved moments and X‘^ + Y^ has the 
same relaxation parameter Sg as the fourth order moment. It is thus consistent for 
this draws to be independent of a. 

We now focus on the moments ® : the parameter a appears only in the third order 
moments. For the TRTi, X^ + W and the third order polynomials are relaxed with 
the same relaxation parameter Sg. Choosing 

XY^ + a(X^ + y^), X^Y + a(X^ + Y^), 
is then equivalent to choose 

XY^,X^Y, 

and the scheme does not depend on a as the left draw of the figure shows it. For 
the TRT2, X^-hY^ and the third order moments are relaxed with different relaxation 
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Figure 4. Draw of as a function of a for the scheme relative 
to u = V with the moments ([^. Left: TRTi, Sg = 2 — 2 “"* and 
s^ = 2- 2-^. Right: TRTa, Sg = 2 - 2"”^ and 8^ = 2- 2-”. 


parameters: it is expected to have a dependence on a, excepted for the BGK case 
{m = n = 7) involving only one relaxation parameter. 

We now do the same job for the scheme relative to u = V. The figure |^is associated 
with the moments ([^ and the figure]^ with the moments (|^. 

The stability of the scheme relative to u = V depends on a whatever the moments 
(figure Q. The maximum is reached for a = 0 whatever the choice of s. For the 
moments ([^, the stability of the TRTi is not linked to a (figure]^ on the left side). 
Instead, this parameter is influential for the TRT2 (figure]^ on the right side): a = 0 
still corresponds to the optimum. 

We interpret the figure [^corresponding to the moments (j^. Because = \^X and 
= X^Y on the velocity set, relaxing the relative moments associated with (j^ is 
equivalent to relax 

1, X, y, - y2, xy, P^{u, a), a), F|(S, a), 

where 


Pq{u, a) = XY^ + a{-3u^X^ + (A^ - 3(S^)2)X + u^{X^ - 
¥j{u, a) = X^Y + ai-Su^Y^ + (A^ - 3(S^)2)y + u^X^ - (S^)^)), 


(14) 

(15) 
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Figure 5. Draw of as a function of a for the scheme relative 
to u = V with the moments (|^. Left: TRTi, Sg = 2 — 2 “"* and 
s^ = 2- 2-^. Right: TRTa, Sg = 2 - 2"”^ and 8^ = 2- 2-”. 


and 


(16) a) = + 1^ ^ ^ 6(S2^)2)y2 

+ 2u^{-X^ + A{u^f)X + 2W{-\^ + A{uyf)Y 

- 3(S^)2(A2 - {u^f) - 3(S^)2(A2 - {uyf)) . 


Let’s observe the equivalent class of the third order moment given by (14). The 


dependence on a of the stability comes from the term —3au^X . Indeed, relaxing 
(14) is equivalent to relax XY^ — 3avYX‘^ since the moments corresponding to the 


polynomials 1 and X are conserved by the collision. On the contrary, the moment 
associated with is not conserved. For the TRTi, it is a linear combination of the 
moments + Y^ and X^ — Y^ associated with different relaxation parameters Sg 
and s^. For the TRT2, it is associated with Sg whereas ig corresponds to 8p. 

These remarks justify the introduction of the moments (1^ to study the influence of 
the non conserved components X^ and Y^. The figur e [sj implying the moments ([^ 
gives similar results as its analogous for S = 0 (figure pF: the same interpretation is 
still valid. Note that for a = 0, the areas are bigger with u = V (figure]^ than with 
u = 0 (figure]^. This confirms the observations of the section [tT 































14 


FRANgOIS DUBOIS, TONY FEVRIER, AND BENJAMIN GRAILLE 


3. Stability for the Kelvin-Helmholtz instability 

The purpose of this section is to conhrm on a non linear test case the previous linear 
stability results: this test case is the Kelvin-Helmhotz instability [?, ?]. 


We compare six versions of the relative velocity D2Q9 scheme to study the influence 
of the moments, of the velocity held u and of the equilibrium. We consider the 
scheme associated with a = 0 relative to S = 0 and u = u (the huid velocity) for the 
equilibria ([^ and ( |10[ ). We compare it to the choice a = 1 for the relative velocities 
0 and u with the equilibrium (|^. We choose not to consider the product equilibrium 
(10) for a = 1, this equilibrium being introduced for the moments of the cascaded 
scheme [?]. We work with the TRTi dehned by (11): unless otherwise specihed, Sg 
et Sjj are hxed by 

A^AfJe 

5 —, i' = —^, 


where Ug = 1/Sg — 1/2 and cjj^ = l/s^, — 1/2, so that the viscosities /r and u are set to 
0.0366 and lO”"^. 


We test the stability of the scheme by increasing the velocity U dehning the initial 
shear layers 


u^{x,y,0) 


Utanh{k{y - \)) if y ^ ^ 
U tanh(A:(| - y)) if y > ^ 


{x,y) E [0,1]^ 


u^{x, y, 0) = U6 sin(27r(x + i)), (x, y) E [0,1]^. 

This velocity U is chosen as Ma/\/3 for MaE M the Mach number. The parameters k 
and 5 controlling the width of the shear layers and the magnitude of the initial data 
are set to 80 and 0.05. 


We hrst validate the vorticity draws obtained in [?, ?, ?] using the scheme relative 
to the fluid velocity u for the second order truncated equilibrium ([^. This vorticity 
is dehned by 

UJ = d^Uy — dyUx- 

For this simulation, the domain is constituted of 128 x 128 points, the Mach number 
is hxed at 0.04 (A is chosen as in [?] so that U = 1). The hgures and are the 
vorticity plots at time t = 0.6 and t = 1. 

We now present a stability analysis depending on the different parameters for A = 1. 
We expect to conhrm the linear stability results. The scheme is considered stable 
if it has not broken after 2000 iterations. The table |6] contains the maximal stable 
Mach number Ma for different meshes at 0.01 close. The table [^presents the greater 
Reynolds number Re = 1/u stable at 1000 close for different meshes and Ma= 0.09. 
Since we discuss on the Reynolds number, the viscosity u becomes a parameter. 







0.2 0.4 0.6 0.8 


Figure 7. Vorticity draw at a t=l. 


We obtain results consistent with the linear stability study. First, choosing a scheme 
relative to u = u has a positive effect if a = 0, negative if a = 1. We must choose the 
moments of the D 2 Q 9 cascaded scheme to improve the stability. This improvement 
occurs whatever the equilibrium and the mesh: the stability limit = 2 is stable 
(table and high Mach numbers are reached for this scheme (table |^. Second, the 
d’Humieres scheme is independent of a as for the linear stability study. Its stability 
area is smaller than the scheme relative to u = u when a = 0 , greater when a = 1. 
Finally, the equilibrium does not influence a lot the stability unlike the linear case. 
The obtained values are close whatever the choice of the equilibrium. It is impor¬ 
tant to note that the tablej^exhibits a convergence of all the schemes as Ax decreases. 
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Space step Ax 

1/16 

1/32 

1/64 

1/128 

1/256 

1/512 

1/1024 

Corresponding Sg 

0.44 

0.25 

0.13 

0.07 

0.03 

0.02 

0.01 

Corresponding 

1.98 

1.96 

1.93 

1.86 

1.73 

1.53 

1.24 


0.18 

0.13 

0.12 

0.12 

0.12 

0.12 

0.12 

0.96 

0.82 

0.62 

0.49 

0.43 

0.39 

0.39 

0.18 

0.13 

0.12 

0.12 

0.12 

0.12 

0.12 

0.92 

0.80 

0.62 

0.50 

0.43 

0.39 

0.39 

0.18 

0.13 

0.12 

0.12 

0.12 

0.12 

0.12 

0.09 

0.07 

0.06 

0.05 

0.05 

0.05 

0.05 


a = 0, M = 0, equilibrium (9 I 
a = 0, u = u, equilibrium (91 


a = 0, M = 0 , equilibrium (101 


a = 0, u = u, equilibrium (10) 


a = 1, M = 0, equilibrium (9 I 
a = 1, u = u, equilibrium (91 


Table 6. Maximum of Ma stable according to the mesh. The two 
last columns indicates a convergence of each scheme as Ax decreases. 


Space step Ax 

1/16 

1/32 

1/64 

1/128 

Corresponding Sg 

0.44 

0.25 

0.13 

0.07 

Corresponding 

1.98 

1.96 

1.93 

1.86 


a = 0, M = 0 , equilibrium 
a = 0, u = u, equilibrium 
a = 0, M = 0 , equilibrium 
a = 0, u = u, equilibrium 
a = 1, u = 0, equilibrium 
a = 1, u = u, equilibrium 


(91 

(91 


(lOl 


(91 

(91 


s^ = 2 

21.10^ 

17.10'^ 

17.10^ 

s^ = 2 
s. = 2 

s^=2 

21.10^ 

s^ = 2 
17.10^ 

s^ = 2 
17.10^ 

s^ = 2 
s^ = 2 

^ = 2 
21.10^ 

s^ = 2 
17.10^ 

s^ = 2 
17.10^ 

10.10^ 

6.10^ 

4.10^ 

4.10^ 


Table 7. Maximum of the Reynolds number stable for Ma = 0.09 
according to the mesh. 


We now characterize the behaviour of the scheme when the diffusion is weak (when 
the relaxation parameters are close to 2). The table presents the maximal Ma 
stable for decreasing bulk viscosity /r. The domain is constituted of 128^ points and 
u is still equal to 10~^. 


This table is also consistent with the linear stability study. When Sg and are far 
from each other, the linear case (the tables and presents an important gain for 
the scheme relative io u = u for a = 0 whatever the equilibrium. These results are 
confirmed by the first columns of the table [^corresponding to take a big bulk viscos¬ 
ity. When ^ tends to 0, the stability areas for the different choices of u are expected 
to be close at fixed equilibrium: this case corresponds to close parameters Sg and 
regime where the linear stability results are homogeneous in u. This behaviour is 
confirmed by the table [^ indeed, the four cases associated with the equilibrium (|^ 
have the same stability areas when the bulk viscosity is smaller than 10“^. Similarly, 
the two cases associated with the equilibrium (10) have close stability areas for these 
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Viscosity /r 

10"^ 

0.5.10-2 

10"^ 

10-^ 

10-® 

Corresponding Sg 

0.23 

0.41 

1.13 

1.86 

1.98 

Corresponding s^, 

1.86 

1.86 

1.86 

1.86 

1.86 


a = 0, M = 0, equilibrium 
a = 0, u = u, equilibrium 
a = 0, u = 0, equilibrium 
a = 0, u = u, equilibrium 
a = 1, u = 0, equilibrium 
a = 1, u = u, equilibrium 


(91 

(91 

m 


(lOl 


(91 

(9) 


0.22 

0.29 

0.43 

0.38 

0.32 

0.73 

0.76 

0.68 

0.63 

0.60 

0.22 

0.30 

0.45 

0.38 

0.32 

0.72 

0.76 

0.76 

0.63 

0.61 

0.22 

0.29 

0.43 

0.38 

0.32 

0.10 

0.14 

0.32 

0.38 

0.32 


Table 8. Maximum of Ma stable according to fj,. 


u 


0 = 0, equilibrium 
0 = 1, equilibrium 
0 = 0, equilibrium 


(91 

(91 

(lOl 


0 

0 . 2 m 

0.4 m 

0 . 6 m 

0 . 8 m 

U 

1 . 2 m 

1.4 m 

0.12 

0.15 

0.21 

0.34 

0.60 

0.49 

0.42 

0.33 

0.12 

0.11 

0.09 

0.07 

0.06 

0.05 

0.05 

0.04 

0.12 

0.15 

0.21 

0.34 

0.60 

0.50 

0.42 

0.33 


Table 9. Maximum of Ma stable according to u. 


viscosities. 

The table deals with the influence of the velocity held u on the stability: other 
choices than 0 and u are considered. We determine the maximal Mach number stable 
for different u according to the choice of the moments. We study the two choices 
0 = 0 and a = 1 for a mesh of 128^ points. 

This table is an evidence of the importance of the moments for the relative velocity 
schemes. Taking a velocity different from 0 provides stability improvements only 
for a = 0. These moments improve the numerical stability for u = u compared to 
the d’Humieres scheme whatever the equilibrium. Instead, choosing u ^ 0 for the 
moments ([^ deteriorates the stability of the scheme. The most stable choice for 
a = 1 corresponds to the d’Humieres scheme. 

4. Conclusion 

We have studied the numerical stability of the relative velocity D2Q9 scheme with two 
conservation laws. A linear stability study was presented and strenghtened by a non 
linear test case for the compressible Navier-Stokes equations: the Kelvin-Helmholtz 
instability. The main conclusion of the article is the following: the relative velocity 
schemes improve or deteriorate the stability of the d’Humieres schemes and it depends 
strongly on the choice of the moments. An improvement occurs if the moments of 
the cascaded scheme are chosen whatever the equilibrium. It is bigger when one 
viscosity is very small and the other is important. The usual set of moments and its 
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orthogonalized analogous deteriorates the stability of the d’Humieres scheme. This 
degradation originates from the presence of second order components in the third 
and fourth order moments. These components do not appear for the moments of the 
cascaded scheme that explains the better stability behaviour. 
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